DIAGRAMA CASUAL DEL CASO
library("deSolve")
#La variable exógena growth population rate se calculó utlizando la siguiente fórmula:
#Anual population rate= 1-(final value/starting value)^1/n
#Donde:
#final value = 2
#starting value = 1
#n = 408 años
#1-(2^(1/408))-1 = 0.0017
parameters<-c(
growth.population.rate=0.0017, #[% per year]
emigration.ratio=0.05, #[% per year]
consumed.food.per.person=400, #[kg per person per year]
intensity.of.agriculture=1,
x=1.3595 #Valor calculado, ver NOTA 1.
)
InitialConditions <- c(
population=100000, #[persons]
fertility.of.agricultural.land=5000000, #[kg/km2]
agricultural.land=8, #[km2]
forest.land=5000 #[km2]
)
#-1000: El signo negativo representa el tiempo B.C, por lo cual es el año 1000 B.C.
#2000: El signo posistivo representa el tiempo A.D, por lo cual es el año 2000 A.D.
times <- seq(-1000,2000,1)
#Definimos método de integración
intg.method<-c("rk4")
#Especificamos modelo
MAYAS <- function(t, state, parameters) {
with(as.list(c(state,parameters)), {
#auxiliary variables
food.produced <- agricultural.land*fertility.of.agricultural.land #[kg]
demand.food <- population*consumed.food.per.person #[kg]
food.gap <- demand.food - food.produced #[kg]
people.without.food <- food.gap/consumed.food.per.person #[persons]
#flow variables
natural.population.increase <- growth.population.rate * population #[persons]
emigration <- emigration.ratio*people.without.food #[persons]
fertility.losses <- fertility.of.agricultural.land* min(2,(agricultural.land/forest.land)^x)/intensity.of.agriculture #[kg/km2]
deforestation <-min(food.gap/max(fertility.of.agricultural.land,1),forest.land/4)/intensity.of.agriculture #[km2]
#state variables
dpopulation <- natural.population.increase - emigration #[persons]
dfertility.of.agricultural.land <- (-1)*fertility.losses #[kg/km2]
dagricultural.land <- deforestation #[km2]
dforest.land <- (-1)*deforestation #[km2]
list(c(dpopulation, dfertility.of.agricultural.land, dagricultural.land, dforest.land))
})
}
out.normal <- ode(y = InitialConditions,
times = times,
func = MAYAS,
parms = parameters,
method =intg.method)
datos <- data.frame(out.normal)
NOTA 1 ¿Cómo se estimo la variable x?
library(ggplot2)
#Gráfica de población
population.collapse <- ggplot (datos, aes(x = time, y = population)) +
geom_line(color="skyblue3")+
xlab('Time (years)') +
ylab('Mayan Population (persons)') +
ggtitle('Collapse of Mayan population') +
theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
population.collapse
#Gráfica de fertility.of.agricultural.land
ggplot (datos, aes(x = time, y = fertility.of.agricultural.land)) +
geom_line(color="darkgoldenrod2")+
xlab('Time (years)') +
ylab('Fertility of agricultural land (kg/km2)') +
ggtitle('Behaviour: Fertility of Agricultural Land') +
theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
#Gráfica de agricultural.land
ggplot (datos, aes(x = time, y = agricultural.land)) +
geom_line(color="chocolate1")+
xlab('Time (years)') +
ylab('Agricultural land (kg/km2*year)') +
ggtitle('Behaviour: Agricultural Land') +
theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
#Gráfica de forest
ggplot (datos, aes(x = time, y = forest.land)) +
geom_line(color="darkgreen")+
xlab('Time (years)') +
ylab('Forest Land (km2)') +
ggtitle('Behaviour: Forest Land') +
theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
library("deSolve")
library(ggplot2)
#Normal
parameters1<-c(
growth.population.rate=0.0017,
emigration.ratio=0.05,
consumed.food.per.person=400,
intensity.of.agriculture=1,
x=1.3595
)
#Cambio marginal 1
parameters2<-c(
growth.population.rate=0.0017,
emigration.ratio=0.08,
consumed.food.per.person=400,
intensity.of.agriculture=1,
x=1.3595)
#Cambio marginal 3
parameters3<-c(
growth.population.rate=0.0017,
emigration.ratio=0.1,
consumed.food.per.person=400,
intensity.of.agriculture=1,
x=1.3595)
InitialConditions <- c(
population=100000,
fertility.of.agricultural.land=5000000,
agricultural.land=8,
forest.land=5000 )
times <- seq(-1000,2000,1)
# Definimos método de integración
intg.method<-c("rk4")
#Especificamos modelo
MAYAS <- function(t, state, parameters) {
with(as.list(c(state,parameters)), {
#auxiliary variables
food.produced <- agricultural.land*fertility.of.agricultural.land
demand.food <- population*consumed.food.per.person
food.gap <- demand.food - food.produced
people.without.food <- food.gap/consumed.food.per.person
#flow variables
natural.population.increase <- growth.population.rate * population
emigration <- emigration.ratio*people.without.food
fertility.losses <- fertility.of.agricultural.land* min(2,(agricultural.land/forest.land)^x)/intensity.of.agriculture
deforestation <-min(food.gap/max(fertility.of.agricultural.land,1),forest.land/4)/intensity.of.agriculture
#state variables
dpopulation <- natural.population.increase - emigration
dfertility.of.agricultural.land <- (-1)*fertility.losses
dagricultural.land <- deforestation
dforest.land <- (-1)*deforestation
list(c(dpopulation, dfertility.of.agricultural.land, dagricultural.land, dforest.land),emigration.ratio=emigration.ratio)
})
}
out1 <- ode(y = InitialConditions,
times = times,
func = MAYAS,
parms = parameters1,
method =intg.method)
out2 <- ode(y = InitialConditions,
times = times,
func = MAYAS,
parms = parameters2,
method =intg.method)
out3 <- ode(y = InitialConditions,
times = times,
func = MAYAS,
parms = parameters3,
method =intg.method)
datos1 <- data.frame(out1)
datos2 <- data.frame(out2)
datos3 <- data.frame(out3)
#Gráfica de población
ggplot() +
geom_line(data = datos1, aes(time,population),
color = "darkgreen")+
geom_line(data = datos2, aes(time,population),
color = "orange")+
geom_line(data = datos3, aes(time,population),
color = "blue")+
labs(x = "Time (years)", y = "Mayan population (persons)")+
ggtitle("Mayan population") + theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
#Gráfica de fertiliy of agricultural land
ggplot() +
geom_line(data = datos1, aes(time,fertility.of.agricultural.land),
color = "darkgreen")+
geom_line(data = datos2, aes(time,fertility.of.agricultural.land),
color = "orange")+
geom_line(data = datos3, aes(time,fertility.of.agricultural.land),
color = "blue")+
labs(x = "Time (years)", y = "Fertility of agricultural land (kg/km2)")+
ggtitle("Behaviour of fertility of agricultural land (kg/km2)") + theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
#Gráfica de agricultural land
ggplot() +
geom_line(data = datos1, aes(time,agricultural.land),
color = "darkgreen")+
geom_line(data = datos2, aes(time,agricultural.land),
color = "orange")+
geom_line(data = datos3, aes(time,agricultural.land),
color = "blue")+
labs(x = "Time (years)", y = " Agricultural land (km2)")+
ggtitle("Behaviour of agricultural land (km2)") + theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
#Gráfica de forest land
ggplot() +
geom_line(data = datos1, aes(time,forest.land),
color = "darkgreen")+
geom_line(data = datos2, aes(time,forest.land),
color = "orange")+
geom_line(data = datos3, aes(time,forest.land),
color = "blue")+
labs(x = "Time (years)", y = " Forest land (km2)")+
ggtitle("Behaviour of forest land (km2)") + theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
Descripción de la gráfica: El color verde de la gráfica representa el comportamiento de las variables de estado del sistema con una growth population rate del 5%, el color naranja representa el comportamiento con una growth population rate del 8%, mientras que la línea de color azul representa el comportamiento con un 10% de growth population rate. Con base a las gráficas podemos decir que el sistema no es sensible a cambios marginales en los párametros.
Cambio marginal en el parámetro consumed.food.per.person
library("deSolve")
#Normal
parameters1<-c(
growth.population.rate=0.0017,
emigration.ratio=0.05,
consumed.food.per.person=400,
intensity.of.agriculture=1,
x=1.3595
)
#Cambio 1
parameters2<-c(
growth.population.rate=0.0017,
emigration.ratio=0.08,
consumed.food.per.person=400,
intensity.of.agriculture=0.8,
x=1.3595)
#Cambio 2
parameters3<-c(
growth.population.rate=0.0017,
emigration.ratio=0.1,
consumed.food.per.person=400,
intensity.of.agriculture=1.5,
x=1.3595)
InitialConditions <- c(
population=100000,
fertility.of.agricultural.land=5000000,
agricultural.land=8,
forest.land=5000
)
times <- seq(-1000,2000,1)
# Definimos método de integración
intg.method<-c("rk4")
#Especificamos modelo
MAYAS <- function(t, state, parameters) {
with(as.list(c(state,parameters)), {
#auxiliary variables
food.produced <- agricultural.land*fertility.of.agricultural.land
demand.food <- population*consumed.food.per.person
food.gap <- demand.food - food.produced
people.without.food <- food.gap/consumed.food.per.person
#flow variables
natural.population.increase <- growth.population.rate * population
emigration <- emigration.ratio*people.without.food
fertility.losses <- fertility.of.agricultural.land* min(2,(agricultural.land/forest.land)^x)/intensity.of.agriculture
deforestation <-min(food.gap/max(fertility.of.agricultural.land,1),forest.land/4)/intensity.of.agriculture
#state variables
dpopulation <- natural.population.increase - emigration
dfertility.of.agricultural.land <- (-1)*fertility.losses
dagricultural.land <- deforestation
dforest.land <- (-1)*deforestation
list(c(dpopulation, dfertility.of.agricultural.land, dagricultural.land, dforest.land))
})
}
out1 <- ode(y = InitialConditions,
times = times,
func = MAYAS,
parms = parameters1,
method =intg.method)
out2 <- ode(y = InitialConditions,
times = times,
func = MAYAS,
parms = parameters2,
method =intg.method)
out3 <- ode(y = InitialConditions,
times = times,
func = MAYAS,
parms = parameters3,
method =intg.method)
datos1 <- data.frame(out1)
datos2 <- data.frame(out2)
datos3 <- data.frame(out3)
library(ggplot2)
#Gráfica de población
ggplot() +
geom_line(data = datos1, aes(time,population),
color = "dark green")+
geom_line(data = datos2, aes(time,population),
color = "orange")+
geom_line(data = datos3, aes(time,population),
color = "blue")+
labs(x = "Time (years)", y = "Mayan population (persons)")+
ggtitle("Mayan population") + theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
#Gráfica de fertiliy of agricultural land
ggplot() +
geom_line(data = datos1, aes(time,fertility.of.agricultural.land),
color = "darkgreen")+
geom_line(data = datos2, aes(time,fertility.of.agricultural.land),
color = "orange")+
geom_line(data = datos3, aes(time,fertility.of.agricultural.land),
color = "blue")+
labs(x = "Time (years)", y = "Fertility of agricultural land (kg/km2)")+
ggtitle("Behaviour of fertility of agricultural land (kg/km2)") + theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
#Gráfica de agricultural land
ggplot() +
geom_line(data = datos1, aes(time,agricultural.land),
color = "darkgreen")+
geom_line(data = datos2, aes(time,agricultural.land),
color = "orange")+
geom_line(data = datos3, aes(time,agricultural.land),
color = "blue")+
labs(x = "Time (years)", y = " Agricultural land (km2)")+
ggtitle("Behaviour of agricultural land (km2)") + theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
#Gráfica de forest land
ggplot() +
geom_line(data = datos1, aes(time,forest.land),
color = "darkgreen")+
geom_line(data = datos2, aes(time,forest.land),
color = "orange")+
geom_line(data = datos3, aes(time,forest.land),
color = "blue")+
labs(x = "Time (years)", y = " Forest land (km2)")+
ggtitle("Behaviour of forest land (km2)") + theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
Descripción de la gráfica: El color verde de la gráfica representa el comportamiento de las variables de estado del sistema con un consume food per person 400 kg por persona al año, el color naranja representa el comportamiento con un consume food per person 200 kg por persona al año, mientras que la línea de color azul representa el comportamiento con 600 kg por persona al año, de consume food per person. Con base a las gráficas podemos decir que el sistema no es sensible a cambios marginales en los párametros.
Para evitar el colapso de la población se decició implementar la política pública denominada “Fertilization workshop”, la cual impacta directamente a la variable de estado fertily of agricultural land lo cual permite que la demand of food no crezca más que la food produced y por tanto la población maya no tenga que emigrar. Todo lo anterior se puede observar en el CLD presentado a conitnuación, así como el comportamiento de la población (en color morado) que no muestra un colpaso a diferencia del comportamiento representado en color azul, donde la population colapsa en el año 800 A.D.
#With political
parameters<-c(
growth.population.rate=0.0017,
emigration.ratio=0.05,
consumed.food.per.person=400,
intensity.of.agriculture=1,
x=1.3595,
fertilization.per.km= 1 #km2
)
InitialConditions <- c(
#Siempre debe serguir el mismo orden cuando se definen.
population=100000,
fertility.of.agricultural.land=5000000,
agricultural.land=8,
forest.land=5000
)
times <- seq(-1000,2000,1)
# Definimos método de integración
intg.method<-c("rk4")
#Especificamos modelo
MAYAS <- function(t, state, parameters) {
with(as.list(c(state,parameters)), {
#auxiliary variables
food.produced <- agricultural.land*fertility.of.agricultural.land
demand.food <- population*consumed.food.per.person
food.gap <- demand.food - food.produced
people.without.food <- food.gap/consumed.food.per.person
#flow variables
natural.population.increase <- growth.population.rate * population
emigration <- emigration.ratio*people.without.food
fertility.losses <- fertility.of.agricultural.land* min(2,(agricultural.land/forest.land)^x)/intensity.of.agriculture
deforestation <-min(food.gap/max(fertility.of.agricultural.land,1),forest.land/4)/intensity.of.agriculture
fertilizers <- population*fertilization.per.km
#state variables
dpopulation <- natural.population.increase - emigration
dfertility.of.agricultural.land <- fertilizers- fertility.losses
dagricultural.land <- deforestation
dforest.land <- (-1)*deforestation
list(c(dpopulation, dfertility.of.agricultural.land, dagricultural.land, dforest.land))
})
}
out.political <- ode(y = InitialConditions,
times = times,
func = MAYAS,
parms = parameters,
method =intg.method)
data.political <- data.frame(out.political)
#Gráfica de población
population.political <- ggplot (data.political, aes(x = time, y = population)) +
geom_line(color="purple")+
xlab('Time (years)') +
ylab('Mayan Population (persons)') +
ggtitle('Behaviuor population with Policy') +
theme_minimal()+
theme(plot.title = element_text(hjust = 0.5))
library(ggpubr)
comparation <- ggarrange(population.collapse, population.political,
ncol=2,nrow=1
)
comparation